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ABSTRACT 


Numerical  evaluation  of  the  optimum  estimate 
via  configural  sampling  involves  evaluation  of 
several  double  integrals.  These  integrals 
represent  expectations  over  a  distribution  condi¬ 
tioned  on  the  observed  configuration.  Theoreti¬ 
cally,  any  location-scale  invariant  definition  of 
the  configuration  will  suffice,  though  numeri¬ 
cally,  some  choices  are  hotter  than  others.  A 
related  concern  is  the  change-of-var iables  used  to 
map  the  region  of  integration,  originally  the 
half-plane,  onto  a  fixed  region,  such  as  the  unit 
square.  This  report  is  of  use  to  the  reader  both 
as  a  guide  to  the  pitfalls  and  curiosities  of  the 
computations  presently  recommended,  and  as  an 
addendum  to  Technical  Reports  185  and  190  [see 
references]  on  the  configural  polysampling 


1.  Introduction. 

As  discussed  in  Technical  Report  187  (Pregibon  and 
Tukey,  1981),  configural  polysampling  techniques  are  useful 
in:  (1)  determining  the  maximum  attainable  efficiency  in  a 

particular  sampling  situation,  (2)  determining  the  maximum 
attainable  polyefficiency  in  a  particular  polysituation,  and 
(3)  guiding  the  modification  of  a  robust  estimate  with  the 
aim  of  increasing  its  polyefficiency.  In  section  2,  we 
describe  the  procedure  and  computations  involved  in  (1) 
above.  Section  3  discusses  those  involved  in  (2).  Item  (3) 
requires  assessing  the  behavior  of  an  estimate  at  particular 
data  configurations  and  will  not  be  discussed  in  this 
report.  An  appendix  lists  the  programs,  including  the  FOR¬ 
TRAN  integrator,  used  in  the  computations. 

2.  Single  situations. 

*  background  * 

Consider  a  sample  { x . : i=l , . . . ,n}  from  a  particular 
situation  {f^  i=l,...,n}  where  the  f.  are  location-scale 
densities.  Following  Bruce,  Pregibon  and  Tukey  (1981),  the 
situation  is  termed  simple  if  fj  =  f  for  all  i;  otherwise  the 
situation  is  termed  compound.  For  example 
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is  a  simple  situation  whereas 


n-1  X's  -  Gau(0,l) 
one  X  ~  Gau (0 , 100 ) 
is  a  compound  situation. 

Configural  methods  require  transformation  from  the 
observed  sample  to  its  location-scale  invariant  representa¬ 
tion.  In  most  cases,  the  configuration  is  expressed  via 

transformation  of  the  order  statistics  y,  <  y„  <•••<  y  . 

l  —  z  —  —  n 

The  general  form  of  the  chang e-o f-var iables  is 
r  =  r(y) 
s  =  S(_£) 

ci  =  i  =  l , . . .  ,n  , 

where  r  is  a  measure  of  location  and  s  a  measure  of  scale. 

Configural  methods  restrict  attention  to  location-scale 
invariant  estimators  t (^)  =  t(r+S£)  =  r+st(c) .  This  allows 
the  determination  of  the  minimum  mean  squared  error  (MSS) 
estimate  of  location  conditional  on  the. observed  configura¬ 
tion  { c .  :  i»l ,  . . .  ,n}  .  Without  loss  of  generality,  assume 
that  fj  is  centered  at  *j=0  with  scale  cr=l.  Then  the  condi¬ 
tional  mean  squared  error  of  the  estimate  is 
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MSE  { t  (y)  I  c}  =  E  {  r+  st  (  c)  I  c }  ^  . 

This  quantity  is  minimized  by 
tQ(c)  =  -E  {rs  !  c  }/E  (  s2  |_c  } 

with 

MSE{tQ(_y)  |  c)  =  E{rs|c}  tQ(c)  +  E{r2|c}  . 

Averaging  MSE{tQ(y)l£)  over  the  distribution  of  configura¬ 
tions  provides  an  estimate  of  the  uncond itional  variance  at 
tQ(£).  The  estimate  tQ(^;)  is  uncond  i  tional  ly  minimum  vari¬ 
ance  for  a  symmetric  situation  since  in  that  case  the  uncon¬ 
ditional  bias  is  zero.  For  any  particular  sample,  the 
optimal  estimate  and  its  conditional  MSE  can  be  computed  by 
numerical  evaluation  of  the  conditional  expectations  as  we 
now  describe. 


*  computational  details  * 

Samples  are  generated  in  a  subroutine,  and  passed  in 
common  to  the  main  program.  The  program  is  shown  in  listing 

I  in  the  appendix.  The  data  may  correspond  to  a  sample  from 
either  a  simple  or  compound  situation.  A  sorting  subroutine 
(listing  4  in  the  appendix)  provides  the  order  statistics 

y,  < . .  .<  y  . 

I I  —  —  1  n 

The  configuration  { c ^ }  is  formed  by  making  the  change 
of  variables 
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-  4  - 


r  =  r(v) 


s  =  s(v) 


4  =  (yi-r)/s 


.n-2 


The  Jacobian  of  this  transformation  is  s  Thus,  in  terms 

of  our  new  coordinates,  we  have  the  probability  element 

f(y)dy  =  sn  2f { r+sc) drdsdc  . 

The  marginal  density  of  £  is 


h(c)  =  J  f  sn  2f(r+sc)drds  . 

r  s 


The  range  of  integration  in  this  expression  is  the  half¬ 
plane.  In  order  to  improve  the  accuracy  of  a  fixed-point 
quadrature,  we  map  the  half-plane  onto  the  unit-square  vie 
(see  Relies  and  Rogers,  1977): 

l 

u  =  l/(l+exp(n2 ( log  s-log  s* ) 1 )  0  £  u  £  1 


v  =  l/(l+exp[n  (c-c*') /s]  ) 


0  <  v  <  1 


where  s*  and  r*  are  appropriate  centering  values  for  the 
bivariate  conditional  density 

g(r,s|c)  =  sn-2f { r+sc) /h ( c)  . 

The  Jacobian  of  this  transformation  is 
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. -  - - --nny  |-Ii  ar-  -  |]  I  irttfUMnititfriifi-ir  If-  1 


J 


5 


„  ,  ,  s*  ,  1-u  4 \  I  n 

J  (u  (v)  =  — g-(-jj~) 

nu  v 


<rk> 


Thus,  in  terras  of  our  new  coordinates,  we  have  the  probabil¬ 
ity  element 

f(^)dy  =  J(u,v)s(u,v)n  2  f  (  r  (  u ,  v)  +s  (  u ,  v)  c)  dudvdc 


=  g  ( u  ,v  ,£)  dudvd£  . 


*  cubature  * 

The  evaluation  of  the  required  conditional  expectations 
can  now  be  carried  out  by  two-dimensional  numerical  integra¬ 
tion  (cubature).  The  following  integrals  (each  defined  on 
the  unit  square) : 

(1)  h ( c)  =  g ( u,v ,£) dudv 

(2)  E(s2|£)h(£)  =  s(u,v) 2g ( u ,v ,£) dudv 

(3)  E(r2|£)h(£)  =  r ( u ,v) 2g ( u ,v ,£) dudv 

(4)  E ( r s I £) h (£)  =  r ( u  ,v) s ( u , v) g ( u ,v ,£) dudv . 

have  so  far  been  done  using  a  24  point  Gaussian  quadrature 
rule  in  both  dimensions  (but  see  below).  Thus,  for  example, 
(1)  is  computed  as 

24  24 

h(c)  =  55  w.w.  g(z.,z.,c) 

j=l  k=l  3  K  3  K  “ 

where  { w. : i=l , . . . , 24 }  and  { z . : i=l , . . . , 24 }  are  optimally 
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chosen  ’..eights  and  evaluation  points  along  one  dimension. 

In  particular,  these  values  are  chosen  so  that  the  finite 
sum  is  exactly  h(c)  for  one  dimensional  polynomials  g(z)  un 
to  degree  47  ((Krylov,  19^2,  pp. 110-111  and  337-340), 
(Abramowitz  and  Stegun,  1970)).  The  two  dimensional 
integrator  is  exact  for  a  function  g(zl,z2)  such  that 
g(z2|zl),  the  function  conditioned  on  the  value  of  zl,  is  a 
47  degree  polynomial  and  such  thac  the  one  dimensional 
integrals  are  a  47  degree  polynomial.  The  two  dimensional 
integrator  is  thus  exact  for  a  function  g(zl,z2)  which  is 
not  above  degree  47  in  either  of  the  two  variables. 

A  listing  of  the  one  dimensional  Gaussian  quadrature 

subroutine  used  in  the  calculations  is  given  in  the  Aopend.ix 

(listing  5).  Figure  1  shows  the  grid  of  points  (z.,z.)  on 

J  K 

the  unit  square  at  v.'hich  the  bivariate  function  is  evaluated. 

in  the  integration.  Figure  2  shows  the  grid  of  quadrature 

coefficients,  w^w^,  used  in  the  24  point  quadrature.  The 

values  shown  are  the  quadrature  coefficients  for  the  14^! 

points  in  the  quarter-square  (0<z.<.5,  0<z,<.5),  where  each 

3  K 

weight  has  been  multiplied  by  10^.  Note  that  the  weights 
have  been  plotted  on  an  equally  spaced  grid  but  that  the 
weights  shown  in  figure  2  are  associated  with  points  on  the 
unequally  spaced  grid  (figure  1).  (The  numbers  of  the 
points  (1-24)  with  which  the  weights  are  associated  are 
labelled  in  the  figure.)  The  coefficients  for  the  57m 
points  on  the  unit  square  are  derived  from  the  values  shown 
in  figure  2  using  the  fact  that 
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quadrature  coefficient  (.f+c)  =  quadrature  coefficient  (,5-c) 

for  the  values  of  c  used  in  the  quadrature  program.  Figure 
3  shows  the  values  of  loa^0  (w^w^)  +  5 .  These  are  again  plot¬ 
ted  on  an  equally  spaced  grid  but  are  associated  with  the 
points  of  figure  1. 

As  noted,  the  24  point  Gaussian  quadrature  integrates 
polynomials  up  to  47  exactly,  and  was  useful  for  testing 
purposes.  We  anticioate  reducing  the  number  of  points 
evaluated  to  fewer  than  575  for  post-testing  computations. 

The  two  dimensional  integration  is  obtained  by  provid¬ 
ing  the  integrator  a  function  which  is  itself  a  one¬ 
dimensional  integral.  In  essence,  the  subroutine  calls 
itself.  However,  since  recursive  function  calls  are  not 
supported  in  Fortran,  the  subroutine  must  invoke  a.  copy  of. 
itself  compiled  under  a  different  name.  The  function  argu¬ 
ment  of  the  call  to  the  copy  of  the  integrator  does  the 

actual  functional  evaluations  g(z.,z.  ,c).  A.s  each  of  the 

3  k  — 

integrals  ( 1 ) —  ( 4 )  has  kernel  g(u,v,£),  the  24x24  grid  of 
values  of  g(z^,z,<r£)  need  only  be  computed  once.  We  take 
advantage  of  this  property  by  storing  the  matrix  g(z.,zk,£) 
after  evaluation  of  (1),  and  using  these  values  for  evalua¬ 
tion  of  (2)  -  (4).  This  provides  us  with  the  quantities 

h (£)  =  (1) 

E(S2|£)  =  (2 ) / ( 1 ) 
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-  - - - —^t****** 


E(r2|c)  =  (3  )  /  ( 1  > 

E ( rs| c)  =  (4 ) / ( 1 ) 

as  are  needed  in  calculating  tQ(_c)  and  VSE  ( t  (y)  I  c)  . 


The 

output  from  a  typical 

run  of  th 

e  prog  ran 

is  a 

( M  +  l )  x  7 

h  (£1 ) 

array  of 
E(s2ICj) 

the  form: 

E ( rs 1 ) 

E ( r 2 ! £x ) 

V-Sp 

W 

MSB (t(^j |Cj ) 

. 

. 

. 

. 

• 

. 

. 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

h<£N) 

E(s2|cn) 

E (rs! cM) 

E (r2|cv) 

U\ 

o 

4J 

W 

”SE(t(vN,) lcM) 

h(c) 

E  (s2) 

E  (rs) 

E  (r2) 

V*5 

‘•'SE  (t(v)  ) 

Each  of 

the  first 

N  rows,  corresponds  to 

estimates 

o  f  the 

conditional  expectations  given  an  individual  con f ig ur a t ion . 
The  final  row  provides  the  estimates  of  the  unconditional 
expectation  obtained  as  the  average  over  configurations. 

*  major  choices  * 

There  are  several  choices  in  the  computational  pro¬ 
cedure  outlined  above  which  have  an  effect  on  the  accuracy 
of  the  results.  These  include  the  choice  of  r*  and  s*  an^ 
the  forms  of  r (y)  and  s(^)  used  in  the  transformation  from 
the  data  to  the  configuration.  We  now  discuss  these. 

Relies  and  Rogers  (1977)  use  the  transformation  (r,s) 

-•»  (u,v)  where  r*  and  s*  are  the  points  at  which  the  density 
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_  o  _ 

n-2  . .  . 

s  • f ( r+sc ) 

attains  its  maximum.  They  state  that  this  transformation 
causes  the  functions  we  integrate  to  be  more  closely  con¬ 
stant  on  their  domains.  To  reduce  computation  time  and 
expense,  it  appears  advantageous  to  choose  r*  and  s*  by  a 
method  other  than  that  suggested  by  Relies  and.  Rogers. 


The  possibility  of  using 


has  been  tested  for  various  functional  forms  of  r  and  s. 
The  form 


r  -  y(l) 
s  =  y(n)-y(l)  , 

i.e.,  the  minimum  as  the  location  estimate  and  the  range  of 

the  data  as  the  scale  estimate  has  the  property  of  putting 

the  configuration  on  the  interval  [0,1],  However,  when 

these  forms  are  used  with  r*=r  .  and  s*=s  ,  ,  the  estimates 

ods  obs 

produced  for  some  samples  are  very  inaccurate. 

The  problem  with  this  approach  can  be  seen  in  a  close 
look  at  the  integration  for  a  straggling  sample.  Samples 
with  large  values  of  y(n)-yd)  were  observed  when  the  data 
were  generated  from  the  slash.  The  density  of  the  slash  is 
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in 


—  ^(l-exp{-  iy2})  for  y/0 

\l2iry2  2 

— - —  for  y=n  . 

2  \  |  2  ir 

Alternately,  slash  is  defined  as  the  ratio  of  an  independent 
Gaussian  to  a  uniform  (0,1)  random  variable.  The  slash  den¬ 
sity  is  like  the  Gaussian  in  the  middle  and  like  the  Cauchy 
in  the  tails,  and  so  has  much  longer  tails  than  the  Gaus¬ 
sian. 


In  samples  with  a  large  range  the  contribution  from 
several  points  on  the  24x24  grid  used  in  the  quadrature 
swamp  all  others  and  the  double  integration  reduces  to  the 
weighted  sum  of  the  values  of  the  function  at  only  a  few 
points.  Figure  (4)  shows  the  24x24  grid  of  powers  of  10  1 
of  the  values  of  g(u,v,c)  used,  for  a  particular  configura¬ 
tion,  in  evaluating  the  double  integral.  This  plot  is  for  a 
sample  of  n=20  with  y ( 2 0 ) — y ( 1 )  =  41,722,  and  with  y ( 1 5 ) — y ( 5 ) 
=  2.808.  The  values  here  and  in  the  figures  5-9  are  shown 
on  an  equally  spaced  grid,  but  correspond  to  points  on  the 
grid  shown  in  figure  1. 

As  an  alternative,  the  location  and  scale  measures 

r  *  midpivot  =  mean  of  the  pivots 

s  =  pivotspread  =  difference  of  the  pivots 

were  tried  and  used  with  r*=r  .  and  s*  =  s  .  in  the  second 

obs  obs 

transformation.  The  pivot  depth  is  defined  as  the  integer 
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in 


11 


part  of  the  hinge  depth,  i.e.,  for  sample  size  n,  pivot 
depth  =  w^ere  t^e  brackets  denote  integer  part. 

The  pivots  are  then  the  order  statistics  with  depth  =  pivot 


depth  and  the  midpivot  is  the  average  of  the  two  pivots. 
Figure  (5)  shows  the  24x24  grid  of  powers  of  10  1  of  the 
values  of  the  function  g (u,v ,c) /J (u,v)  (i.e.,  without  the 

Jacobian  J(u,v)  from  the  second  transformation)  when  these 


new  values  are  used.  Figure  (S)  is  the  comparable  plot  for 
the  function  g(u,v,c).  Comparing  figures  (4)  and  (^),  we 
see  a  much  more  constant  order  of  magnitude  of  the  function 
over  the  domain  when  the  midpivot  and  pivotspread  are  used. 
Figure  (7)  shows  the  grid  of  powers  of  10  *  of  the  values  of 
the  product  of  the  function  g(u,v,c)  and  the  quadrature 
coefficients  used  in  the  integration. 


In  an  attempt  to  make  the  surface  we  integrate  over 
still  more  constant,  r*  and  s*  were  moved  to  correspond  to 
(*)  in  figure  (5).  The  results  are  shown  in  figure  (8), 
where  the  values  plotted  on  the  grid  are  again  powers  of 
10  *  for  the  function  values.  Also  of  interest  is  the 
change  in  the  optimum  estimate  for  the  original 


r*  -  r  .  =  midpivot 

obs 

s*  =  sol3s  »  pivotspread 

and  the  relocated  r*  and  s* .  The  estimate  values  are 
-.70761(38  and  -.  7074012,  respectively.  This  small  change 
(.0002156)  in  the  values  of  the  estimate  leads  us  to  ques- 
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tion  the  gain  fron  recentering. 

3.  Bisampling. 

In  the  previous  section,  computations  for  the  case  when 
data  are  generated  from  and  used  as  if  they  are  from  the 
same  situation  were  described.  In  bisampling  (see  listing  2 
in  the  appendix),  we  distinguish  between  the  generating 
situation  and  the  evaluating  situation.  The  former  is  the 
situation  actually  generating  the  data,  while  the  latter  is 
the  situation  we  treat  the  data  as  being  from  and  at  which 
we  evaluate  the  optimum  estimate. 

Suppose  we  have  two  situations,  for  example,  slash  and 
Gaussian,  f  and  f  _ .  We  generate  a  sample  from  the  Gaussian 

S  VJ 

and  proceed  as  described  in  the  previous  section  to  calcu¬ 
late  the  minimum  variance  estimate  for  the  associated  confi¬ 
guration.  Here  the  Gaussian  is  both  the  generating  and  the 
evaluating  distribution.  We  then  use  the  same  data  and  con¬ 
figuration  and  treat  it  as  being  generated  by  slash,  i.e.  we 
have  a  Gaussian  generating  and  a  slash  evaluating  situation. 
A  similar  procedure  is  followed  with  generated  slash  data. 

In  bisampling,  we  also  calculate  weights,  wG  and  wg  as 

wG  *  fG  (c) /(c(G*  fG  (c) +dg*f  s  (c)  ) 

w s  =  f s ( £) /(dG • f G (c) +ds  *  f s ( c) ) 

where  <*_  and  <X  are  the  sampling  fractions,  N_/(N^+N  )  an<* 
ij  s  j  y  s 

N  / (N_+N  ),  for  the  Gaussian  and  slash,  respectively,  w  is 
s  u  s  o 
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the  weight  pr opor t ional  to  the  probability  that  the  confi¬ 
guration  is  Gaussian  given  that  the  configuration  is  one  of 

N,.  Gaussian  configurations  or  M  slash  configurations;  w  is 
-j  s  s 

defined  similarly  for  the  slash.  These  weights  are  used  in 
calculating  the  average  MSE^  and  MEEs. 

The  output  from  this  program  is 
"G  Eg{s2^'  t0(i’>  *SEG 

ws  Es(s2l£l  t“(i)  WSEs 

for  each  of  the  samples  from  the  Gaussian  and  for  each  of 
the  samples  from  the  slash.  Evaluation  of  the  maximum 
attainable  biefficiency  (for  slash  and  Gaussian  data)  using 
this  output  is  presently  under  consideration  (see  listing  3 
in  the  appendix  and  (Tukey,  1981a)). 

4.  Conclusions. 

Computing  the  optimum  estimate  for  a  situation  using 
configural  sampling  or  configural  polysampling  methods 
involves  the  evaluation  of  several  double  integrals.  The 
choices  of  (1)  functional  forms  of  r  and  s  used  in 
transforming  the  data  to  the  configurations,  and  (2)  the 
values  of  r*  and  s*  as  appropriate  central  values  of  the 
bivariate  density  f  (r,s)  affect  the  precision  and  accuracy 
of  the  numerical  integrations.  The  choices 
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r  =  nidpivot 


s  =  pivotspread 


and 


r*  =  r 


obs 


s*  =  s 


obs 


give  well-balanced  functions  and,  thereby,  good  integral 
evaluations  and  estimates.  This  choice  also  keeps  computa¬ 
tion  costs  to  a  reasonable  level  and  below  those  of  some 
alternative  choices. 
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Figure  8:  Powers  of  10  for  the  values  of  g(u,v,£)  at 

which  the  function  is  evaluated  for  integration 
(r*  and  s*  relocated). 
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GTVT'*- *TE”F  ( Tlir)  *TE,,PP  ( jwr )  **J *Trvpc  ( tvp)  ** v 

PETUm 

FMF 


nonnnnnnno 


A  5 


LTSTIFC  2 

CALCULATE?  T«E  NECESFAPY  COUPLE  INTEGRALS,  OPTIMAL 
ESTIMATE  T{Y)  ,  NSE {T (Y) )  ,  AND  WEIGHTS,  FIP<=T  F^F 
GAUSSIAN  CONFIGURATIONS  EVALUATES  r ?  GAUSSIAN  'NT 
SLASH  CONFIGURATION? ,  THEN  FOR  SLASH  CCNFTGUP' TICKS 
EVALUATED  AS  GAUSSIAN  A  AT  SLASH  CONFTCUP  ATTCM*7 . 

CAN  EE  EXPANDED  TO  CNC  ,  ICCISTIC,  AND  CAUCHY. 

CTPER  COMMENTS  AS  It,’  L  IS  TNG  1. 


IMPLICIT  FE^L*?  (  A.-H ,  C-2 ) 

REAL*?  Cd”)  ,.AA’S(12) 

COTTON  /APE a  ] /C,N/AFFA  4/PFT, SST/APEA 2/L  ,F  ,  J / a F r *  3/V/A PE A.5/INC 
gy'Tt'NAL  FTNT,FINTn 
csnrc=i327i 1 . r'1 
PEArts,!??)  N,vs*.’-p 
15"  FCP?'AT  (IT?, IT"*) 

pr.«;  a'*  psr p  j  £*?  ( NC  A.N P) 

EC  SEE  LL=1 ,3,2 
PC  LOOP=l  ,N?’.NP 

c^ll  r.Ai7rEV{LL.,rsrrr) 

C*LL  SOFT?  (C ,?•’) 
rr:TE(P,43T)  (C(i),t«i,v> 

NT-’*  (N+l )  /? 

N'T  1  =  (NN+1  )  /2 

NN2“N-NNl+] 

f-(C(k::2)+c(piu)  )/2.r? 

C*C(NN2)-C(NN!) 

pj'  2C"  "*1  v 

2S?  C(T)=(C(T)-F)/C 

FST=r 

C*t  C  ^  -5  C 

?U;>?.F3 
EC  35r  L=1 , 3 , 2 

ccns-ccnst (r ) 

IF  (L.  EQ. 1  )  INC*? 

IP(L.E2.3)  INC-4 

K=f 

J*" 

CALL  CTNT2A ( F7NT, FSC ) 

K  =  2 
J-'1 
ZKr** 

call  DTNT24  ( FINT'* ,  F2C ) 

F2n=F20/F'?r 
K  *1 

Jal 

INF*" 

CALL  DTNT24 { , FI  I ) 

FI  3  *F2  1  /?'*n 
K*" 

J*2 


- 


A  9 


35'’ 


i  /in'! 

41'* 

45" 

c;r» 


c 


1C 

15 

2C 

25 


in 

ic 


Tfinsfl 

CaLL  DINT2* (FIMTC , Fp?) 

FC2=FC2/F?" 

F""  =FPC  *COFS 
PTTC  =  -F11  /F2*'* 

PTTX*P+S*PTTC 
PIT”  'r  F=  F 1 1  *PITC+F:’2 
A?’0-  (irC4-!  )  ■  Fi*  1 
? !’”  (INC-*  2)  *F2'2 
A?’?(TtsC+3)*P:mX 
A  vc  ( ifrc+4  )  *PITM?E 

CF'sftji'  +  Ffli 

CCt’TIf’UE 

"T.’F  (1)  *AK?  (1  )/FV'< 

?vc  ( 5 )  =AT5  ( 5 )  /cVr 

I’FTTF  (P  ,  4'’" )  ( A  vc  (Jr)  ,17  =  1  <  p ) 

FTP” A" (3 (4  Tlc . 7/) ) 

FCF”*T  PEI*.  7) 

^Q..mT»rrfC 
Q  <r*  rn  i  »7  y  jr 

5T0P 

rvj 

rCUPI.F  PPECT9I0K  FITVCTTOM  FTFT(V) 
TrPITCIT  P r  *  P  (A-P,C-") 

EVTEF”?  L  CT”'7 , C7>;,Tf 
CC””CH  /A°FA3/VV 
W=V 

C  a  1 1  A.™??*  (GItf7,AtT5) 

FI!'7*A.K3 

PETUPK 

EFTFY  FT’TTn(V) 

W=V 

C  A I  £  ATNT24  (GT11T'’ ,  A  t’S) 

FTt’T(?*AvS 

PFTL’PP 

fmt 

FUsPCVTTV.t  F3?’rEV(LL ,  OFFED 
TMPI ICTT  PEA I *P  ( A -P , C-" ) 

TEA  I  *9  C<lf*'1) 

CC’T’O!!  /APFA  !/C  ,P 

P»TR  PI/3.1415C2FP?57P28/ 

GT  TO  (lf",20,?',,*a,5^)  , TT. 
rc  15  i*i 
C ( I ) *GGt’QF ( ESEFr) 

PE"'UPN 

P'*  25  T*1  ,V 

C  ( I)  *GGV?F  (rcFED 

C(»4)-70.rn*C(F) 

PETUFF 
rc  35  1*1, N 

c ( i ) -cgvqf (crrrr) /cgheff (dsfed 


£r 

45 

50 

ce 


2 

•5 


1C 

15 

2C 


tc 


3n 


pcrupy 

CO  ^  T  =  1 ,V 
?  cr=jrjinpi  (C'-EFC) 

C(I)*r'L',G(’»FG/(l.r'’,-PpC)  ) 

PETUP.K 
ro  55  1-3  ,13 

C(I)*CTAK  (PI*  (GC-U5FF  (D5EEC)  -0 .50'*)  ! 

EETURK 

EFC 

rCL’-iE  PFrcrsiCF  fuectiof  coreT(u 

IMPI ICTT  PE.*L*°  (O-HjC-Z) 
pr & L * 0  C(1 O'*) 

co:::*oi:  /»ffm/c  ,r 

r s  t:t  ^3 . 1  *  1  c92F5"3  c7(52£>/ 

CC!,rT*l .  ro. 

GO  TC  (1,1, 1,3,2)  ,T. 

pri^/jNC ''*£>•"  (  o  _  r'"1  * p  j  ) 

^QTicrf<— cti*  *v 

cor:^T*(l.ro/Pi)  **r 

PETUPM 

ryir 

pprjpy.p  ppnGIcTOM  Ft ~''.r'VTQ™  CT,Irn(F) 

IMPLICIT  PE’L*?  (5-r,C-o) 

FEM.*0  z(]<*'')  ■J’Pf’PS  ( 57c  )  ,TE’*?e(57D  ,mEVF(c75) 

cc"r:oK/?.rr*  i/?  f  /?cm,c'rc/i!'pri2/t ,  f  ,  j  /pee*.  -,/v/f  pe*  5  /TFr 

TVpsTVPi] 

E,''arFEO?  ,T'  (D 
P0T’*1 .  rfl/ri?CpT  (OF) 

■yr'-pe  (T»:p.)  *c«pT*  (  (  (1  .  /D  **PCD 

TEF.PP  (IFF)  *TEMPr  (T?T)  *POr*ELCG(  (1  .rn-v)  /V)  +Fr,T' 

TC'P  (Tt’F)  *DN*CtCG (Tf  'PS  (IFF) )  -ELOC  (CF*C*V*  (1  .CC-rM  *  (1  .Ff-V) ) 

ru’-'1  =c .  r° 

grr»‘  2»o  ep 

GO  TOP  ,20 ,30 ,4 '1,5!’)  ,1 
CO  15-  T«i  ,M 

X=TP'PF  ( 7!’r )  +TE’*Pc;  ( !VE )  *  Z  ( I ) 

FL’0 1  *OUC 3 -y*X/2  .  co 
GO1  TC  ^ A 
EC  25  1*3 ,K 

X*TE!'?p  (IFF)  +TE** PE  ( IFF)  *Z  ( I ) 

X2E2*X*X/2.rp 
Evpc^*X2r  2*3  .  f'onr' 

3F(FypCv  .  GT .  I”.-. CD  EXPC»?*3  7C  .r« 

Fl"'l*cUMl-X2E2 

EU”2*Sr”.2+EEXP  (Tv PCD  /(FV'*1F  .  rD 
£UV1*SCM1+CLCC  (rU’!2) 

GO  TO  6? 

CO  35  1*1, V 

Y •TC" PP ( IFD )  +TE”pr'  (P-T)  *2(1) 

x2=x*x 


.  ....  — 


A8 


EXPCN=2  .  5En*X.2 

IF(EXPOK  .CT.  173. En)  EXPCr=]  1*.  DO 

sr:2*i .  co-dexp (-exp cr) 

G U”  ] =T  UN 1 +C 1 CG ( CU" 2 /X  2 ) 

GC  TC  SO 
DC  45  1=1  ,i: 

X=TENPF.  ( INC )  +TE*'PS  ( I  IT'D )  *2(1) 

expci:=x 

IF  (CAES  (EXPO?:)  .Cl.  170  .C")  EXPCV=F«  ICN  ( 1 72  .  r*  ,  EXPON) 
rUV2  =  l.rP+EEX?  (EXPCl: ) 


45 

SU”1*SU** 

GC  TC  63 

50 

EC  55  1= 

X=TE'*PP( 

55 

£E!'.  1  =CL"’1 

60 

T'"P=TC”  P 

IF(D?FS  (7?-:P)  .CT.  17".r(’)  TCP-mCN  (1 7* .  rn  ,  TUP) 
TE.”P  ( I!JE )  =DEXP  (TKP) 

GINT*TE'-'P(I!-TD) 

FETUFN 

ENTRY  GINTO (U) 

TND=TND+1 

C-INTO=TEUP  ( IND)  *TE*'PP  ( T!'E )  **J*TEMT>5  ( INF)  **!' 
FCR'!5T  (CIS .  “) 

FETUP" 

cvr. 


i 


C  tI5TT?T  3 

C  C7*  ICUL*TEC  p TEFFIC Ir^fn  (rLT  C!I  f  £c,>mT,,•  f 

c  sn*rcv  prices,  **ir  excefs  (srr  ^ukfY/IQ0!^ ' 

C  FCP  IflR  S^MPir?  E*CH  op  CM/SSI.7**1  ^wr  SLAStT  (**.*Mprr  ?n). 

IMPLICIT  RE.7'  I  *°  (  *  — H  ,  o-Z ) 
riMENSIO”  a  (2)  ,F-(2)  ,V(2) 

F.V2*  EPF ,  ITL/1 .  C-"-*  ,  lr/ 

CO  fcn  K*1 ,2 

**i  _  pf 

vc  =  ".rr 

FVG  =  ^.C° 

pvo  *  P.rd 
cc 

C  pr*  ^  0/"?F'rOr’P? ■•"Tp**  VPTCrf'Fp  ayp  pT'T'*»^y  VBpr*nrr? 

fc^c ( ^ ,  11 )  vc,Fu:,vr,rid 

li  For:  *•?(///,  c?  5.  ?2y,  n  ci5. 7,  ??v,r.i  p.7) 

FVG  *  PVG  +  VG*FHl 
PVE  =  PV?  +  WS* 3112 
VI  =  VI  +  vc 
V2  *  V2  +  V? 
ff£  CCI’TZVIT 

C  FVC  (PV0)  12  TFr  "EIC!?TEn  ’‘VEPJ'GF  CF  pT"'i,>”  VA ria’,Cr'?  FCF 

C  GAUEEIA'.T  ( «•  t  a  “ li )  . 

TF(K  . cr .  1)  C COO  F2n 
FVG  =  (PVG /VI  +  FO/I.r'1 
PVS  *  (PVG/V2  +  PF )  /2  •  P'1 
CPTC  F5'7 

62C  PG  *  PVC/V1 
PC  *  PVS/V2 
f-50  ccrTiruc 

rr.ITE(F,2?)  PVG,  PVS 
22  FCF”5T(4rlf . 7) 

23  F0F2  *7  (1  °X  ,  '  THE  PIT:*PK  V?  FI2.VCEF  I FE :  ,,7CI5.'7) 

FEVIVC  7 
KF*  *  P 
5(1)  =  2.7 
E(l)  *  2.3 
4  02  CCKCirL’E 
prr-rvr  "* 

RET'I'-T  ? 
r-<">  1  r>  ff  t  —  i 

w‘^  x.’  f  * 

C  C?t GUI* TEE  CP7IKAI  T  »VC  FELACTVE  pypr^c  v^PT^VCCS  FOP  SPFOTFTEr 

C  SH^CCV  FFTCCr  (5  vt  p)  . 

rr?C(7,2)  X5  ,X!F,  i:c,  20,00,  VP,  SF,T2 
1  FOp!!rT  (F3X,riF.7,//15X,EIF.’7,/rl5.7f2r1F.7  ,/nlS  . 7 , 2CI  r  .  7) 

XI  «  (X16  -  X E ) * ( X 1 F  -  X?) 

VK1  *  SG/PVC 
W?:2  «  SS/PVS 

?  *  TG*A.  (1)  *VC*V«T  +  TS*P  ( 1 )  *VS*V?I2 
T  ■  T/(VC*A(1)  *V»H  +  VS*P  ( 1  )  *V*»?) 

AFFO  -  (TG-T)*(TG-T)*VF']/X1 

s  (03~r')  *  (""S-T)  *  VP 7 /X F 
VriTE (2 ,2)  T,VC, Ar«G,VS, APES 


MO 


2  FOPV  AT  ( 5 Dl  6  .  ° ) 

ICO  CONTINUE 
REWIND  8 
CO  25?  K*1 , 2 
VC  *  0. 

VS  *  C. 

UH1  *0. 

V.H2  *  M 

C  CALCULATE?  VETCr^Er  fVEF4CF  CF  RELATIVE  EVCE?0  V*'rT?NCFr  F^r  T 

pn  ’■  —  1  ^  *  (*,**• 

P FAD(P , 2)  T,t:0,  AM?C,NT,  'MSS 
VC  =  VG  +  rr**MSC 

VS  *  V5  +  t;e*ivec 

vhi  =r:-:?  +  wo 
ri-2  =  vr?  +  r? 

202  CONTINUE 

IF  (K  .EO.  1)  GOTO  22  S 

vc  *  (vg/vhi  +vvc ) /? . 

VS  *  (VS/M’2  +  WS)/2. 

COTC  25r 

22"  VVC  *  VC/,‘t,i 
VV?  =  VS/l'H.? 

25^  CONTINUE 
KH  *  KHM 
V(l)  *  VS  -  VO 

white  (MS)  vc, vs 

5  FCF.MAT  ( 1  OX  ,  ’  T!T  EXCESS  VaP.Ir-VC£S  .nct:  ‘,2FIC.") 

C  ITEFMEE,  CHANGING  CHACON  PPICCC  L,T,TIE  (VC-VSK.MM 

C  CP  IP 

IF(CAFS(VG-VF)  .IE.  EPS  .CP.  K!T  .CF.  ITL)  GOTO  5nn 
IF(KH  .EC.  1)  GCTC  2^ 

AH  *  Ml)  -  V(!)*(M1)-M?))/(V(1)-VP)  ) 

JF { aH  .CF.  1.)  a"=  1. 

M2)  =  .Ml) 

E  ( 2 )  *  EU) 

A(l)  =  AH 
P  ( 1 )  =  1.  -  ?!’ 

V  ( 2 )  =  V  ( 1 ) 

GCTC  e"? 

30?  A  (2)  =  Ml) 

F ( 2)  =  E(l) 

V  (2 )  =  V  (1 ) 

I F  (VC  .EC.  DM?  XI  (VG,VP) )  E(l)  «VC/VC*P(M 

m i )  =  i  -  rm 

IF  ( VE  .EQ.  DMA  XI  (VC  ,  Vc )  )  Ml)  *VC/Vf,*M1) 

rm  =  1-*  (l ) 

C-C^O  400 
5??  CONTINUE 

*'PITE  (6  ,  ?)  Ml)  ,P(1  ) 

2  FORMAT (1 OX, 'THE  SHATON  PPICFc  aPEr  ’,2C11.4) 

VFITE (5,4)  VC,V? 

t  F0FM’'T(1FX/ '  E^0FoS  V*  D  .-FT:  CPUSS’  cI’SHr  ,,T'1<:.7) 

END 


r 


All 


C 

C 


IISTIN’C-  £ 

Sl'EFCUTTNE  FOPTP  (V  , »J) 
TUPEICTT  PFAI *8  (A-K 
FCA  L*F  V  (N) 


C  SPELT  SOFT  ALCCFTTH"  CaCM  JULY  19?4 

T  —  1 

X  . 

1  I=T+I 

IF  (I  .LE.  K)  C-C  TC  1 

(.< s  r  _  1 

2  "=v/2 

IF ( V  .EC.  P)  FETUFN 
K  =K-" 

EC  4  J*1 , K 
L=J 

5  I F ( L  .1”.  1)  GO  TC  £ 

IF  (V  (L+?1 )  . CE  .  V(L)  )  CC  TC  4 

X=V(L+r) 

V(L+r-)=V(L5 
V  { L )  *  X 
L=L-*' 

X  *■-'  ~f 

i  CCUTIUUE 

GO  TC  2 
ECO 


■‘hrliii  • 


LISTING  5 


SUE  POUT  IKE  AINT2*  (FCT,y) 

DOUBLE  PPFCISTOF  Y^,c,PCT 
D»TA  A/9 . 5U9/ 

C=.49759',6''999?C1<,*°C'' 

Y=.617C61^909993599^r-2* ( FCT ( »+C) +FCT( A-C) ) 
C*. 4973? 4 277995654 75D5 

Y*Y+. 14 2f 5 6°4  3144 6? 9P2D-1  * ( FCT ( A+C) +FCT ( B-C) ) 
C*. 4 69 13 727f ^"l?f 6  3  809 

Y*Y+.  221  387194  r.e?''pQp7r-l*  (  FCT (»+C)  +FCT  (  B-C)  ) 
C=.4432'177635',22P',S?C7 

Y®Y+  .296 4929 2-1 5771 P  39'* D-l*(  FCT{AXC)  -*-FCr"  ( A-C)  ) 
C*.41C,9B?092°PC95146T''' 

Y*Y+.  26 67324 C7n?e4«i5'»t‘-i*(rCT(Jv^)+FCT(.»-C)  ) 
C*.  3 79 '*6  2095 79027*71  POO 

Y=Y+  .  43'195?P'7 7^507 c^^Pf-1*  (  FCT  ( A.-*-C)  +  F''n*  f  B-C)  ) 

C*.  324f'46925q6P4?77or(’ 

Y=Y+.  49°99326',c27e;60d4c<_i*(  pr-?  (  b+C)  +FC'P  ( A-C)  ) 
C*. 27271 3735f?aaiQ77r« 

Y-Y+.53722135',5~«°2P17C-1*(FC'?'(B+C>  +FCT(»-C)  ) 
C*.216p96753«13r,2257rf’ 

Y=Y+.  57‘752P?Bfl26p?2Ff’10-l*(  FCT(  B  +  r)  +FCT(B-0  ) 
C*.  1575  2 13399  49"  ®1  69C<" 

Y*Y+.699352364f 39^1 6°?D-1*( FCT (B+C)+FCT(A-C) ) 
C=.955596?373*p7P15r'-l 

Y=Y+.P291B7281736]£1 3  Pr-1* ( FCT (A+C) +FC" (A-C) ) 
C*.  320294 4643.1  3^291 3C-1 

Y*Y+.  F'»969e9  7€T3  37<:',78r-l*(  FCT  (A  +  C)  +FCT (B-C)  ) 

RETURN 

END 


Figure  8: 


Powers  of  in  tor  the  values  oi  gvu,v,v-/  o.. 
which  the  function  is  evaluated  for  integration 
(r*  and  s*  reloceted). 


